Assessing the effect of urbanization on regional-scale surface water-groundwater interaction and nitrate transport

Identifying regional-scale surface water-groundwater interactions (SGI) is vital for predicting anthropogenic effects on surface water bodies and underlying aquifers. However, large-scale water and nutrient flux studies rely on surface water or groundwater-focused models. This study aims to model the effect of urbanization, which is usually accompanied by high groundwater abstraction and surface water pollution, particularly in the developing world, on a regional-scale SGI and nitrate loading. In the study area, the urban expansion increased by over 3% in the last decade. The integrated SWAT-MODFLOW model, Soil and Water Assessment Tool (SWAT) and Modular Finite-Difference Groundwater Flow (MODFLOW) coupling code, was used to assess SGI. By coupling SWAT-MODFLOW with Reactive Transport in 3-Dimensions, the nutrient loading to the river from point and non-point sources was also modeled. Basin average annual results show that groundwater discharge declined with increasing groundwater abstraction and increased with Land use/Land cover (LULC) changes. Groundwater recharge decreased significantly in the Belge season (February to May), and the river seepage and groundwater discharge decreased correspondingly. High spatiotemporal changes in SGI and nitrate loading were found under the combined LULC and groundwater abstraction scenarios. The water yield decreased by 15%. In a large part of the region, the nitrate loading increased by 17–250%. Seasonally controlled groundwater abstraction and water quality monitoring are essential in this region.


Study area descriptions
The study area is in Awash River Basin, Ethiopia. Awash River flows from central highlands to the northeast through the rift floor and ends in Lake Abhe, and the basin covers over 10% of Ethiopia's area. A significant part of the basin lies in the Great East African Rift Valley. Compared to other river basins in the country, Awash River Basin is the most developed basin 32,33 . The capital city and several small urban centers, urban-centered industries, small towns, and agricultural activities characterize the upper part of the basin. More than 60% of the country's industries are established in the Awash River Basin 31,34 . Addis Ababa, Mojo, and Adama are among the major industrial areas in the region.
Since most of the industries were established following the river course, the fundamental water quality and quantity-related problems emanate from this arrangement 35 . A few industries have treatment plants 34 . For example, several industries in and around Addis Ababa discharge wastewater with limited or without treatment into the Akaki River 36 . Nevertheless, regardless of the anthropogenic pollution and over-abstraction in some parts of the basin, the water supply of the city depends entirely on the Akaki River Basin 37 . City expansion and population growth are escalating the demand; the contribution of groundwater to Addis Ababa's water supply increased over three folds in the last 15 years 31,38,39 . Groundwater abstraction wells are scattered through the entire basin, including in the city, without protected areas 37 . Despite the vital role of the groundwater in the area, hydrogeological studies and data are extremely scarce or lack a systematic database that can easily be accessed 33,40 .
This study focuses on the most urbanized area of the Awash River Basin. The area is vastly developed, and consequently, the water resources' quality and quantity, both surface and groundwater, are affected 35 . Reservoirs, namely Dire, Legedadi, Gefersa in the upstream area, and Aba Samuel downstream, are the main water bodies in the study area (Fig. 1). Gefersa and Dire reservoirs are relatively small. Aba Samuel reservoir served the city as a water supply source and hydropower generation. The rivers draining Addis Ababa transport high concentration pollutants, and the reservoir became a non-functioning swamp. The major rivers draining the region are known as Little Akaki, Kebena, and Big Akaki. Little Akaki River runs from Menagesha Mountain (upstream of Gefersa reservoir) to Aba Samuel. All three rivers drain the city, but Little Akaki crosses the most populated and industrialized part of the capital city.
Previously, the water supply of Addis Ababa derived from springs around the foot of the Entoto Ridge (northern part of the city) and dug wells in the central and southern parts of the city 38,41 . When the demand continued to increase, dams were constructed to use the rivers as a water supply source. Furthermore, exponential population growth and rapid urban expansion have led to further groundwater exploration in the Akaki well-field in the southern part of the city and several more places around the city. The continuous over-abstraction and lack of protection zone made the aquifer around Akaki well-field delicate 38 .
Waste management in Addis Ababa is in its infancy. The city has a poor sewer system and several residential and commercial buildings discharge their waste directly to the nearby river or runoff drainage network. The only municipal wastewater treatment plant is Kaliti Treatment Plant (KTP). The effluent from KTP is about 7500 m 3 /d with a peak value of 10,000 m 3 /d during rainy periods; it is a negligible amount for a city hosting over four million people 6,38 .
The principal rainy season (Kiremt) in the region lasts from mid-June to September, while the dry season (Bega) is from October to January. In the rest of the months (Belge season), the area gets moderate and  47 and U.S. Geological Survey "earthexplorer" websites, respectively, were used (Fig. 2a, d). The GlobeLand30 LULC classification has ten types. However, in the study region only eight LULC classes, namely, Forest, Cultivated land, Grassland, Shrubland, Wetland, Water Bodies, Bareland, and Artificial surfaces ( Table 1). The soil data were acquired from Harmonized World Soil Database v1.2 48 . Eutric Vertisols cover a large part of the area (Fig. 2b). Daily wind speed, temperature, relative humidity, precipitation, and solar radiation from 1980 to 2013 were collected from the National Metrological Agency of Ethiopia.  www.nature.com/scientificreports/ The study area, Akaki River Basin, was subdivided into 37 subbasins and an HRU of about 1363. The topography varies from 3239 m around Entoto Mountain to 1813 m above sea level near the outlet of the basin. About 77% of the area has a slope of 0-15% (Fig. 2c). The HRU was generated using the zero-threshold value of LULC, soil, and slope option.
Groundwater flow model and hydrogeological descriptions. Topographically the study area favors the emergence of springs in the hills and valleys. The elevation descends from north to south, and more of the central, eastern, and southeastern parts are low areas covered with thick Quaternary deposits. Details of geology and hydrogeology data are important for building a representative model. However, like any other part of the country, data scarcity is the principal limiting factor for groundwater-related studies in particular. And therefore, the model developed for this study relied on prior studies 41,49,50 and low resolution or short-duration data collected from the Ministry of Water Resources, Irrigation, and Electricity (EMWRIE) and Addis Ababa Water and Sewerage Authority (AAWSA).
Miocene-Pleistocene volcanic successions dominate the geology in Akaki River Basin, including acidic and intermediate lava flows, basaltic lava flows, and pyroclastic flows forming an interlayered sequence with quaternary faults 39 . The region has highly variable and complex aquifers. The aquifers in southern Addis Ababa are mainly young volcanic rocks of lava flows and tectonic fractures. The thickness and hydraulic conductivity of the unconsolidated sediments govern the subsurface infiltration, overlying weathered and fractured porous volcanic rocks with relatively high infiltration capacity, and favoring circulation and storage of subsurface water. Alluvial, residual, and lacustrine clay deposits dominate and vary with topography and geomorphology. The alluvial deposits, composed of clay with increasing thickness from north to south, are the predominant geological types along the rivers. A thin layer of residual clay soil covers the ridges and hillsides in and around the city. Thick residual clay soils cover high plains around the basin boundary. Black cotton lacustrine deposits cover the lower part of the region. In the area around the Akaki well-field, thick alluvial, lacustrine clay deposits overlay the underlying volcanic rocks interbedded between the volcanic rocks such as scoria, scoriaceous basalt, and basalts. As conceptualized by Ayenew et al. 49 , the regional groundwater flow is from north to south. More details of the geology and hydrogeology of the basin are presented by several researchers, e.g., 41,49 .
The central volcanic hills and mountains, not crossed by faults, bound the river basin; the hydraulic conductivity is highly variable, 0.01-550 m/day 49 . The primary geologic classes in the basin are basalts, scoria, rhyolites, trachytes, ignimbrites, trachybasalts, and tuff of varying ages (Fig. 3a). The scoria and scoriaceous basalts with primary porosity and permeability make up highly productive aquifers in the region. The highly weathered and fractured basalts, ignimbrite, fractured tuff, and pyroclastics are also highly productive aquifers with secondary porosity and permeability. Fractured basalts, sparsely spaced joints and vesicles, ignimbrites, and agglomerates form moderately productive aquifers. The transmissivity is also highly variable, increasing from the highlands to Akaki well-field 41,49 .
Even though the study area is data-scarce, the hydrogeology framework is described in several previous works, e.g., 49,51 . Previously, the groundwater flow modeling was performed by using two vertical layers and 400 × 400 m cell dimensions 49 . However, the thickness of the top layer is thin compared to the bottom layer, it varies from "a www.nature.com/scientificreports/ few meters to 53 m" 49 . In this study, it was found that the data which identify the thickness of the first layer was insufficient to build a reasonable model. And thus, a simplified convertible single-layer model was built. The mode area was discretized into 300 × 300 m horizontal resolution, which results in 18,976 active cells. The top and bottom boundaries were interpolated from the well-completion reports. The most important groundwater recharge-discharge areas (boundary conditions) are groundwater pumping, springs, reservoirs, river and canal networks, and percolation from precipitation. Even though pumping wells are distributed all over the basin, the primary groundwater pumping area is the Akaki well-field. In the area, pumping wells with varying depths up to 604 m below the surface and a pumping rate of 13,193 m 3 /day were identified and incorporated into the model using the Well Package. Details of the reservoir operation are unavailable, and thus, to simplify the data requirements, it was modeled as General Head Boundary. The major rivers were imported from the SWAT model and modeled using the River Package. Canal or small rivers in a few places were modeled using Drain Package (Fig. 3b). The initial groundwater head was interpolated from the static water level record.
Integrated SWAT-MODFLOW model. Both SWAT and MODFLOW models are widely used and well known among hydrologists and hydrogeologists. And therefore, here, the general equation of SWAT and MOD-FLOW models and a brief description of the integrated SWAT-MODFLOW code are given (Eqs. 1 and 2). The general equation of the SWAT model is as follows 42 : where SW o and SW t are the quantities of initial and final soil water (mm/day), respectively, t is the time (days), PCP is the precipitation (mm/day), Q s is the surface runoff (mm/day), ET is the evapotranspiration (mm/day), Perc is the percolation (mm/day), and Q gw is the flow from the aquifer (mm/day). MODFLOW is physically-based algorithm that solves groundwater flow through porous earth material 43,52 . The model is based on the following general partial differential equation 43 : where K xx , K yy , and K zz are principal components of the hydraulic conductivity tensor in the x, y, and z special directions, W is the source or sink, S s is the specific storage (1/L), h is the hydraulic head (L), and t is time.
Among the advantages of using the integrated SWAT-MODFLOW is water yield computation. The basic water yield calculation in the SWAT model is based on the following equation 42 : where WY is water yield (L 3 /T), Q s is surface runoff (L 3 /T), Q gw is groundwater flow to the river (L 3 /T), Q lat is lateral soil flow (L 3 /T), T loss is water loss through the riverbed (L 3 /T). T loss , in Eq. 3, is computed with an assumption that when a channel does not receive groundwater, possibly it loses water through the side and bed and joins bank storage or the deep aquifer and approximated using the following relationships: where K ch is the effective hydraulic conductivity of the channel alluvium (L/T), T Q is the flow travel time (T), P ch is wetted perimeter (L), and L ch is channel length (L).
SWAT-MODFLOW uses the MODFLOW River Package that is formulated based on Darcy's law to compute water loss from the river and groundwater discharge to the river with the following mathematical principle: where Q leak is volumetric water flux between the river and the aquifer (L 3 /T); K bed is hydraulic conductivity of the riverbed (L/T); M is riverbed thickness (L); h ch and h gw are river stage and groundwater head (L), respectively; P ch is wetted perimeter (L); L ch is channel length (L). Equation 5 shows that if the groundwater head is higher than the water level in the channel, Q leak will be negative, which means water flows from the aquifer to the channel, and Q leak will be positive if the reverse is true. In SWAT-MODFLOW water yield computation is based on Q leak (Eq. 5) instead of T loss (Eq. 4).
The recent SWAT-MODFLOW algorithm 13,27,44 was used to integrate the calibrated SWAT-2012 and Newton formulated MODFLOW-2005 MODFLOW-NWT: 45 . The current version of SWAT-MODFLOW 44 has been applied to study several hydrological processes, including groundwater recharge 26 , regional-scale SGI 44,53,54 , and the effect of natural and anthropogenic factors on water balance [55][56][57][58][59][60] . However, the diversity of application regions and study objectives is limited compared to the individual SWAT and MODFLOW models.

Reactive transport in 3-dimensions (RT3D) model. Recently, Wei et al. 13 applied SWAT-MODFLOW
for the assessment of non-point source contaminant transport in agricultural areas by integrating it with Reactive Transport in 3-Dimensions RT3D;61 . RT3D is a contaminant and solute transport model that can simulate dispersion, advection, and chemical reactions in a saturated groundwater system 61 . The algorithm has the functionality of solving multi-species reactive transport. The following governing equation describes the fate and transport of aqueous and solid-phase species in multi-dimension saturate porous media 61 : RT3D uses the groundwater hydraulic head, cell-by-cell flow data, and water sources and sinks simulated by the MODFLOW model to establish the groundwater flow field. The detailed linkage processes with SWAT-MODFLOW are provided in 13,16,61 . SWAT-MODFLOW-RT3D combines SWAT-MODFLOW code as a base and RT3D as a sub-routine into a single executable algorithm 13 . The integrated model uses the SWAT 'in-stream' algorithm to rout nitrate through the stream network. In SWAT-MODFLOW, the SWAT model computed recharge is mapped to MODFLOW grid cells, and MODFLOW simulates the groundwater flux in each cell 61 . RT3D uses MODFLOW simulated groundwater flux. RT3D also receives the nitrate concentration from the SWAT model recharge water and river network and computes the change in each grid cell and river mass loading.
In the study area, the major point sources of nitrates are industrial discharges into the river and leaching from waste disposal sites 39,62 . Sources with high concentrations such as industrial discharge and leaching from solid waste disposal sites 63 , along Little Akaki River, were selected (Fig. 3b). In these sites as identified by several researchers, e.g., Angello et al. 63 , the nitrate concentration varies from 33 mg/l to 300 gm/L. The biggest waste disposal site called Koshe was modeled as constant concertation (300 mg/L) point source. The denitrification is specified using single-Monod expression so that the reaction rate depends on the nitrate presence 64 . Thus, the values of the first order-rate denitrification constant (1/d) and Monod half-saturation constant (M/ L 3 ) are required for the model setup. However, these values are not easy to measure and are rarely available for watershed-scale studies, usually assumed based on commonly used values from the literature. In this study area, related studies are unavailable, and therefore, we have referred to worldwide literature, e.g., 15,16,[64][65][66] , and considered 0.05/day and 10 g/m 3 for the first-order denitrification and Monod half-saturation constant. The surface and subsurface hydrological processes, models, and scenario simulations were conceptualized by modifying the schematic diagram from previous studies 13,16,44 . Different colors are used to identify the processes and the corresponding simulating algorithm (Fig. 4).
Model calibration and validation. The SWAT model was calibrated and validated using sequential uncertainty fitting version 2 (SUFI-2) in the SWAT-CUP 67,68 user interface. Monthly river flow at two gauging stations, Little Akaki and Big Akaki, obtained from EMWRIE was used. The calibration and validation performance was evaluated using the Nash-Sutcliffe Efficiency coefficient (NSE) and minimizing the Percent Bias (PBIAS), Root Mean Square Error (RMSE), and coefficient of determination (R 2 ). www.nature.com/scientificreports/ The groundwater model was also calibrated using the Parameter ESTimation (PEST) Package. However, the automatic calibration technique was limited to the steady-state simulation because of the transient groundwater level data scarcity. The important calibrated parameters were hydraulic conductivity, river conductance, and drain conductance. The zonal and pilot point techniques 69 were used to calibrate the hydraulic conductivity using a range value of 0.012 to 375 m/day. The initial groundwater head is a key input for the simulation of integrated surface water and groundwater flow models. In this study, the initial head was based on the calibrated steadystate MODFLOW model simulation. Nevertheless, transient simulation parameters (specific storage, specific yield) of the aquifer calibrated after integration using a trial-and-error technique based on the groundwater level fluctuation, groundwater discharge (or comparing the measured and simulated river flow), and the overall water balance of the study area. In addition to these soft calibration and validation practices previous studies in the region, e.g., 41,49,70,71 , were important sources of initial ranges of parameters and were also used as a benchmark for general the model outputs analogy.
Automatic calibration functions are unavailable in the recent version of SWAT-MODFLOW. Individually calibrated SWAT and MODFLOW for Akaki River Basin were integrated, and moderate manual calibration was performed ( Table 2). In SWAT-MODFLOW-RT3D mainly longitudinal dispersivity and initial concentration of nitrate in groundwater were adjusted by trial and error. We have used longitudinal dispersivity of 120-180 m and constant porosity of 0.3.

Simulation scenarios. Two important water and nutrient mass balance affecting factors, namely LULC
changes focused on the urban expansion and groundwater pumping for urban water supply, were tested using an integrated modeling approach. The effects were evaluated based on the baseline simulation. Therefore, the simulations were baseline, the effect of groundwater pumping (scenario1), urban expansion in the period of 2000-2010 (scenario 2), urban expansion from 2000 to 2010 and groundwater pumping (scenario 3), urban expansion from 2010 to 2020 (scenario 4), and urban expansion between 2010 and 2020 and groundwater pumping (scenario 5). The LULC data is tabulated in "Groundwater flow model and hydrogeological descriptions" Sections ( Table 1). The current groundwater pumping rate was used as a reference and increased by 15-35% to evaluate the effect (Table 3). Related scenarios were simulated in a previous study 49 to assess the aquifer response using MODFLOW. Under all scenarios formulated and simulated in this study, the regional-scale SGI, groundwater recharge, and nutrient loading were the focus of the analysis.

Results and discussion
In 2000, the urban coverage was 14.52% of the total study area. The LULC change assessment shows that the urbanization from 2000 to 2010 was trivial compared to the recent expansion. The Grassland and Artificial surfaces were the most affected classes. From 2010 to 2020, the Grassland decreased by around 1.7%, while the urban area extended by about 3%. The Forrest and Shrubland coverage showed a slight improvement. The effect of these changes on the spatiotemporal distribution of SGI and nitrate seepage from the river to the aquifer is presented compared to the baseline simulation.   72 ranges from "good" to "very good". The peaks and lows of the flow match well with the corresponding basin average precipitation. However, the performance at low and peak flows was dissimilar in the calibration and validation period. The model overestimated the low flows and underestimated the peak flows. For Little Akaki, the validation performance was less than the calibration, while for the Big Akaki, the validation was better than the calibration (Fig. 5). The SWAT-MODFLOW based river flow showed a slight improvement at low flow compared to the SWAT model. In Little Akaki, the flow slightly increased from 1997 onward compared to the preceding years, and rapid urbanization could be one reason for this.
The calibrated groundwater head varied from 1800 to 2653 m (Fig. 6a). Previous research outputs 49 present a closely related calibrated groundwater head. The simulated and observed head agreed fairly with an R 2 value of 0.98 and RMSE of 7.1 (Fig. 6d). The model underestimated the groundwater level in the altitude areas.
The average groundwater discharge to the river map (Fig. 6b) shows that in the study area, the groundwater discharge to the river occurred almost in all sections of the river. However, the seepage from the river to the  www.nature.com/scientificreports/ aquifer was in limited areas, and the magnitude was significantly high, comparatively. The initial nitrate concentration in groundwater was also evaluated using a limited number of measurements obtained from the EMWRIE (Fig. 6c). The measured nitrate concentration varies from zero to 33 mg/L. The initial simulated concentration varied from zero to 0.9 mg/L. The spatial distribution of high and low measured and simulated concentration matches well with an RMSE value of 6.5, even though there is a difference at peak values. While, in a large part of the area, the measured and simulated data show good agreement, the main reason for the peak value of measured and simulated concentration mismatched is that the points with higher concentration have no spatial relation and are limited in number. The spatial distribution of simulated initial nitrate concentration represented the groundwater quality in the study area. To show the model performance in transient groundwater flow simulation (SWAT-MODFLOW), the groundwater discharge to the river and baseflow, analyzed from observed river flow using automated separation techniques 73 , were compared (Fig. 7). The baseflow and groundwater discharge matched fairly, but SWAT-MODFLOW underestimated the discharge in the rainy period and overestimated in dry periods.
The LULC change effect was substantial on average annual basin surface runoff (surq), recharge (rch), and lateral flow (latq). The change in lateral flow from the baseline was quite significant (~ 72%). The river seepage and groundwater discharge were controlled by both groundwater pumping and LULC changes ( Table 5).

Figure 6.
Map showing the calibrated starting groundwater head (a), the spatial distribution of groundwater discharge to the river (b), the initial nitrate concentration (c) in the Akaki Aquifers prepared using the QSWATMOD2 plugin (https:// swat. tamu. edu/ softw are/ swat-modfl ow/), and scatter plot of measured and simulated groundwater head (d). Surface water-groundwater interactions. The effect of the groundwater pumping and urban expansion on SGI in terms of groundwater recharge, river seepage, and water yield over the simulation period was immediately apparent. After several trial-and-error simulations, the effect of groundwater pumping, which increased by an average of the proposed increase rate from the baseline pumping (25%), was considered. The total groundwater abstraction in the baseline scenario was 99,000 m 3 /day. The graphical comparison shows that the effect of groundwater pumping on the monthly water yield was insignificant (Fig. 8a). However, the average value revealed that groundwater pumping caused a decrease in water yield by up to 3%, the LULC change by 9%, and the combined LULC and groundwater pumping by up to 15%. The urban expansion caused to increase in the groundwater recharge in Kiremt and a decrease in the Bega and Belge seasons (Fig. 8b-d). Principally in the Belge season, the recharge showed a significant difference compared to the baseline simulation (Fig. 8d). The recharge under scenarios 2 and 4 had an insignificant difference in all seasons. These results further support the effect of rapid urbanization that occurred from 2010 to 2020 compared to the preceding decade.
Maximum river seepage and groundwater discharge occurred in Belge and Kiremt, respectively (Fig. 9). The LULC change had a decreasing effect on river seepage (Fig. 9a); the effect was significant in the Belge season. On the other hand, the groundwater discharge increased with LULC change (Fig. 9b, scenarios 2 and 4). The groundwater discharge difference between the seasons was trivial. Figure 10 shows the groundwater head differences from the baseline simulation at selected wells. The fluctuation was significant around Akaki well-field. The spatiotemporal groundwater level variability under all scenarios (difference from the baseline simulation) was random and complex. However, comparatively, the groundwater level decreased under almost all simulation scenarios, while in Big Akaki, predominantly in the upstream area, the level showed a slight increase. The central reason is the distribution of the groundwater pumping wells. The results revealed that the response of the aquifers to the LULC and pumping from upstream to downstream is  Table 5. Principal annual basin water balance components, including precipitation (prec), surface runoff (surq), lateral flow (latq), groundwater discharge to the river (gwq), surface water discharge to the aquifer (swgw), and recharge (rch) of the study area under all simulation scenarios (mm). www.nature.com/scientificreports/ highly variable. The groundwater head at the upstream area was more sensitive to LULC changes than groundwater abstraction. The most important flow parameter, which affects groundwater level fluctuation and other water balance components, is groundwater recharge. In addition to the results presented in the previous section, the spatial variability of the groundwater recharge related to the LULC change was analyzed. In the region, the groundwater recharge reached a maximum of 162 m 3 /day, and it varied largely with the LULC change (Fig. 11). Surprisingly, as  . Seasonal surface water discharge to the aquifer (swgw) and groundwater discharge to the river (gwq) differences under the simulation scenarios from the baseline scenario (scenarios-bassline). www.nature.com/scientificreports/ urban expansion increased in the period from 2000 to 2010, the area, which gets moderate to high groundwater recharge, mainly around the boundary of the city, enlarged (compare Fig. 11a, b). The principal reason could be the slight increase in Forest and Shrubland coverage in that period. Compared to scenario 2, the average recharge under scenario 4, particularly in the urban area, increased slightly (Fig. 11c).
Nitrate mass seeping into the aquifer. The effect of urbanization on the average spatiotemporal distribution of nitrate mass seepage from the river to the aquifer is presented in Figs. 12 and 13. The seasonal distribution shows that scenarios 4 and 5 increased the nitrate seepage significantly (Fig. 12a). The seepage was high in the Belge under all simulation scenarios comparatively. These results reasonably mirror the finding on water seepage from the river to the aquifer. However, the peaks and lows trend varied from year to year over the simulation period ( Fig. 12b-d). The potential source for this random nutrient seepage over the years could be the groundwater recharge change corresponding to the LULC changes. The average spatial distribution of nitrate mass seepage, evaluated on the SWAT model subbasin scale, revealed that the seepage was moderate in a large part of the region (600-10,000 kg/day). In most of the areas, the average nitrate seepage increased up to 17-250% from the baseline attributed to the urban expansion. However, in scenario 1, the seepage declined by about 4%. This shows that the average nitrate seepage in the region is more sensitive to LULC change than groundwater pumping. The groundwater pumping effect was significant around Aba Samuel lake, also known as the Akaki well-field (Fig. 13a, b). But the effect had no distinct relationship with the spatial distribution of nitrate seepage in all other scenarios (Fig. 13c-f). A maximum of 116,173 kg/ day mass of nitrate seeps into the aquifer (averaged over the simulation period). Most of the area covering high seepage was in the Little Akaki, the most urbanized subbasin of the region (Fig. 13e, f).

Summary and conclusions
In this work, the effect of urbanization on SGI and nitrate loading was assessed using an integrated SWAT-MODFLOW-RT3D model. The changes in SGI and nitrate seepage related to the LULC change and groundwater pumping were evaluated related to the baseline simulation results. In the simulation of nitrate transport, we have incorporated both point and non-point contaminant sources.
The decadal LULC changes analysis showed that the urban area expanded rapidly in the recent decade. Especially, grassland and artificial surfaces had a continuous change. The water balance was more sensitive to LULC changes compared to groundwater abstraction. Among the principal water balance components, the effect of LULC on later flow was substantial. The annual average water seepage from the river to the aquifer decreased by www.nature.com/scientificreports/ over 4%, while the groundwater discharge decreased by less than 2%. The water yield decreased up to 3% under groundwater abstraction scenarios and up to 9% in LULC change scenarios. Significant changes were found in the Belge season. The combined urban expansion and groundwater pumping compounded the water balance and nitrate seepage changes. The relation between water and nitrate seepage from the river to the aquifer did not exhibit any unique relationship. High nitrate loading was found attributed to LULC changes, and the effect of groundwater abstraction was small. On average, nitrate seepage increased by 17-250% from the baseline in most of the area. This overarching modeling work highlighted the effect of urbanization usually accompanied by an expansion of the built area, groundwater over-abstraction, and surface water pollution, principally in developing countries with poor waste management on groundwater. Even though the focus was not on providing calibrated and validated model, this work relied on limited data, and the results must be interpreted with caution. Incorporating detailed and longtime-recorded data, different contaminant types, and simulation scenarios are suggested for future studies. Further investigation on aquifer recovery and protection strategies is needed for sustainable water resource management in the region.